suppressPackageStartupMessages({
library("here")
library("tidyverse")
library("MplusAutomation")
library("gt")
library("glue")
library("kableExtra")
library("misty")
library("lavaan")
library("AICcmodavg")
library("nonnest2")
library("DiagrammeR")
library("lavaan")
library("tidyLPA")
library("semTools")
library("brms")
library("MBESS")
library("ufs")
library("robmed")
library("careless")
library("psych")
library("BayesFactor")
library("effectsize")
library("tidybayes")
library("emmeans")
library("bayesplot")
library("patchwork")
library("bmlm")
})
options("max.print" = .Machine$integer.max)
# Make random things reproducible
set.seed(1234)
options(mc.cores = 4)
bayesplot_theme_set()
source(here::here("src", "R", "functions", "funs_add_neoffi60_subscales.R"))
source(here::here("src", "R", "functions", "funs_correct_iesr_scores.R"))
source(here::here("src", "R", "functions", "funs_plot_job_qualification.R"))
source(here::here("src", "R", "functions", "funs_generate_all_items_df.R"))
scale_this <- function(x) as.vector(scale(x))
sum_coding <- function(x, lvls = levels(x)) {
# codes the first category with -1
nlvls <- length(lvls)
stopifnot(nlvls > 1)
cont <- diag(nlvls)[, -nlvls, drop = FALSE]
cont[nlvls, ] <- -1
cont <- cont[c(nlvls, 1:(nlvls - 1)), , drop = FALSE]
colnames(cont) <- lvls[-1]
x <- factor(x, levels = lvls)
contrasts(x) <- cont
x
}
all_items <- generate_all_items_df()
There is a problem with IES-R, in the control group. I shift the control distribution of IES-R towards lower values.
temp <- correct_iesr_scores(all_items)
all_items <- temp
ggplot(all_items, aes(x = iesr_ts, colour = is_rescue_worker)) +
geom_density()
all_items |>
group_by(is_rescue_worker) |>
summarize(
avg_iesr = mean(iesr_ts)
)
There are some cases in which the RWs of Toscana and Lombardia did not have the proper qualification. They are assigned to the category of RWs.
all_items$idx <- 1:nrow(all_items)
all_items$is_rescue_worker <- as.character(all_items$is_rescue_worker)
all_items$is_rescue_worker <- ifelse(
all_items$idx < 746, "Si", all_items$is_rescue_worker
)
all_items$idx <- NULL
all_items$commeetee_location <- ifelse(
all_items$is_rescue_worker == "No", "None", all_items$red_cross_commeetee_location
) |>
factor()
all_items$commeetee_location <- ifelse(
all_items$is_rescue_worker == "No" & all_items$age < 25, "students", all_items$commeetee_location
)
all_items$commeetee <- ifelse(
all_items$is_rescue_worker == "No" & all_items$age >= 25, "community_sample",
all_items$commeetee_location
)
all_items$commeetee <-
ifelse(is.na(all_items$commeetee), "community_sample", all_items$commeetee) |>
factor()
sum(is.na(all_items$commeetee))
[1] 0
# all_items <- all_items %>%
# mutate(job_qualification = case_when(
# is_rescue_worker=="Si" & job_qualification == "non_rescue_worker" ~ "team_member",
# TRUE ~ job_qualification))
iesr_ts | trunc(lb = 0) ~ is_rescue_worker + (1 | commeetee),
m0 <- brm(
bf(
iesr_ts ~ is_rescue_worker,
sigma ~ is_rescue_worker
),
family = skew_normal(),
data = all_items,
backend = "cmdstanr"
)
Compiling Stan program...
Warning: incomplete final line found on '/var/folders/hl/dt523djx7_q7xjrthzjpdvc40000gn/T//Rtmp9PpTYj/model-122363e511137.hpp'
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
# m0 <- brm(
# bf(
# iesr_ts ~ is_rescue_worker + (1 | commeetee),
# sigma ~ is_rescue_worker + (1 | commeetee)
# ),
# family = skew_normal(),
# data = all_items,
# backend = "cmdstanr"
# )
pp_check(m0)
Error: object 'm0' not found
summary(m0)
Family: skew_normal
Links: mu = identity; sigma = log; alpha = identity
Formula: iesr_ts ~ is_rescue_worker
sigma ~ is_rescue_worker
Data: all_items (Number of observations: 1068)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 11.72 0.50 10.80 12.78 1.00 1941 2352
sigma_Intercept 2.31 0.04 2.23 2.39 1.00 1966 2109
is_rescue_workerSi 7.16 0.70 5.80 8.54 1.00 1706 1949
sigma_is_rescue_workerSi 0.37 0.05 0.27 0.46 1.00 1828 2083
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
alpha 21.73 2.40 17.27 26.63 1.00 2632 2121
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
me <- conditional_effects(
m0, "is_rescue_worker"
)
plot(me, points = FALSE)
BFt <- BayesFactor::ttestBF(
all_items$ies_ts[all_items$is_rescue_worker == "Si"],
all_items$ies_ts[all_items$is_rescue_worker == "No"],
paired = FALSE
)
effectsize(BFt)
Cohen's d | 95% CI
------------------------
0.69 | [0.55, 0.82]
- Estimated using pooled SD.
Supported families are: ‘acat’, ‘asym_laplace’, ‘bernoulli’, ‘beta’, ‘beta_binomial’, ‘binomial’, ‘categorical’, ‘com_poisson’, ‘cox’, ‘cratio’, ‘cumulative’, ‘custom’, ‘dirichlet’, ‘dirichlet2’, ‘discrete_weibull’, ‘exgaussian’, ‘exponential’, ‘frechet’, ‘gamma’, ‘gaussian’, ‘gen_extreme_value’, ‘geometric’, ‘hurdle_cumulative’, ‘hurdle_gamma’, ‘hurdle_lognormal’, ‘hurdle_negbinomial’, ‘hurdle_poisson’, ‘info’, ‘inverse.gaussian’, ‘logistic_normal’, ‘lognormal’, ‘multinomial’, ‘negbinomial’, ‘negbinomial2’, ‘poisson’, ‘shifted_lognormal’, ‘skew_normal’, ‘sratio’, ‘student’, ‘von_mises’, ‘weibull’, ‘wiener’, ‘zero_inflated_asym_laplace’, ‘zero_inflated_beta’, ‘zero_inflated_beta_binomial’, ‘zero_inflated_binomial’, ‘zero_inflated_negbinomial’, ‘zero_inflated_poisson’, ‘zero_one_inflated_beta’
The sk, ch, mi sub-scales are coded so that high values indicate high self-compassion levels. The sj, is, oi sub-scales are coded so that high values indicate low self-compassion levels.
The ts_sc score has been computed by reversing the coding of the items of the sj, is, oi sub-scales (so that they indicate the absence of self-judgment, absence of isolation, absence of over-identification).
scs_subscales <- with(all_items, data.frame(sk, ch, mi, sj, is, oi, scs_ts))
cor(scs_subscales) |> round(2)
sk ch mi sj is oi scs_ts
sk 1.00 0.52 0.58 -0.39 -0.28 -0.24 0.71
ch 0.52 1.00 0.49 -0.01 -0.03 -0.04 0.45
mi 0.58 0.49 1.00 -0.19 -0.33 -0.35 0.66
sj -0.39 -0.01 -0.19 1.00 0.67 0.66 -0.75
is -0.28 -0.03 -0.33 0.67 1.00 0.80 -0.78
oi -0.24 -0.04 -0.35 0.66 0.80 1.00 -0.78
scs_ts 0.71 0.45 0.66 -0.75 -0.78 -0.78 1.00
In the COPE scale only two factors are identified.
all_items$pos_reinterpretation <- with(all_items, cope_1 + cope_29 + cope_38 + cope_59)
all_items$mental_disengagement <- with(all_items, cope_2 + cope_16 + cope_31 + cope_43)
all_items$venting <- with(all_items, cope_3 + cope_17 + cope_28 + cope_46)
all_items$seeking_instrumental_support <- with(all_items, cope_4 + cope_14 + cope_30 + cope_45)
all_items$active_coping <- with(all_items, cope_5 + cope_25 + cope_47 + cope_58)
all_items$denial <- with(all_items, cope_6 + cope_27 + cope_40 + cope_57)
all_items$religion <- with(all_items, cope_7 + cope_18 + cope_48 + cope_60)
all_items$humor <- with(all_items, cope_8 + cope_20 + cope_36 + cope_50)
all_items$behavioral_disengagement <- with(all_items, cope_9 + cope_24 + cope_37 + cope_51)
all_items$restraint <- with(all_items, cope_10 + cope_22 + cope_41 + cope_49)
all_items$seeking_emotional_support <- with(all_items, cope_11 + cope_23 + cope_34 + cope_52)
all_items$substance_use <- with(all_items, cope_12 + cope_26 + cope_35 + cope_53)
all_items$acceptance <- with(all_items, cope_13 + cope_21 + cope_44 + cope_54)
all_items$suppr_competing_activities <- with(all_items, cope_15 + cope_33 + cope_42 + cope_55)
all_items$planning <- with(all_items, cope_19 + cope_32 + cope_39 + cope_56)
Create COPE sub-scales scores using all items – note that SEM analyses suggest to drop some of the items.
all_items$active_coping <- with(
all_items, pos_reinterpretation + active_coping +
suppr_competing_activities + planning + restraint +
seeking_instrumental_support + acceptance
)
all_items$avoidance_coping <- with(
all_items, mental_disengagement + denial + humor +
behavioral_disengagement + substance_use + religion
)
all_items$soc_emo_coping <- with(
all_items, seeking_instrumental_support +
seeking_emotional_support + venting
)
plot(density(all_items$scs_ts))
pp_check(fit_1)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
me <- conditional_effects(
fit_1, "is_rescue_worker"
)
plot(me, points = FALSE)
summary(fit_1)
Family: student
Links: mu = identity; sigma = log; nu = identity
Formula: scs_ts ~ is_rescue_worker
sigma ~ is_rescue_worker
Data: all_items (Number of observations: 1068)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 74.13 1.01 72.17 76.08 1.00 4836
sigma_Intercept 2.88 0.04 2.81 2.97 1.00 4420
is_rescue_workerSi 7.45 1.16 5.15 9.67 1.00 4717
sigma_is_rescue_workerSi -0.13 0.05 -0.22 -0.03 1.00 5247
Tail_ESS
Intercept 2949
sigma_Intercept 2748
is_rescue_workerSi 2725
sigma_is_rescue_workerSi 2970
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
nu 38.43 15.26 16.30 75.45 1.00 4067 3086
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
BFt <- BayesFactor::ttestBF(
all_items$scs_ts[all_items$is_rescue_worker == "Si"],
all_items$scs_ts[all_items$is_rescue_worker == "No"],
paired = FALSE
)
effectsize(BFt)
Cohen's d | 95% CI
------------------------
0.42 | [0.29, 0.56]
- Estimated using pooled SD.
rw_df <- all_items |>
dplyr::filter(is_rescue_worker == "Si")
rw_df <- rw_df %>%
mutate(job_qualification = case_when(
job_qualification == "non_rescue_worker" ~ "team_member",
TRUE ~ job_qualification))
lpa_scales <- c(
"is_rescue_worker",
"neuroticism", "extraversion", "openness", "agreeableness", "conscientiousness",
"active_coping", "avoidance_coping", "soc_emo_coping",
"iesr_ts",
# "avoiding", "intrusivity", "hyperarousal",
# "sk", "ch", "mi", "sj", "is", "oi",
# "pos_sc",
# "neg_sc",
# "ts_sc",
"mpss_tot"
# "ptgi_total_score"
# "relating_to_others",
# "new_possibilities",
# "personal_strength",
# "appreciation_of_life",
# "spirituality"
)
lpa_rw_df <- subset(rw_df, select=lpa_scales)
# lpa_rw_df <- subset(all_items, select=lpa_scales) |>
# dplyr::filter(is_rescue_worker == "Si")
lpa_rw_df$is_rescue_worker <- NULL
```r
lpa_rw_df %>%
scale() %>%
estimate_profiles(1:6,
variances = c(\equal\, \varying\),
covariances = c(\zero\, \varying\)
#package = \MplusAutomation\
)
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
Model Classes AIC BIC Entropy prob_min prob_max n_min n_max BLRT_p
1 1 21200.56 21292.85 1.00 1.00 1.00 1.00 1.00
1 2 20604.04 20747.10 0.71 0.87 0.94 0.36 0.64 0.00
1 3 20322.46 20516.27 0.75 0.82 0.93 0.18 0.56 0.00
1 4 20180.65 20425.23 0.74 0.83 0.86 0.10 0.42 0.00
1 5 20066.42 20361.76 0.74 0.79 0.88 0.08 0.34 0.00
1 6 20007.30 20353.40 0.73 0.76 0.86 0.06 0.30 0.00
1 7 19941.06 20337.92 0.76 0.77 0.85 0.03 0.31 0.00
1 8 19877.85 20325.48 0.78 0.74 0.88 0.02 0.30 0.00
1 9 19826.21 20324.60 0.78 0.75 0.91 0.02 0.33 0.00
1 10 19793.42 20342.58 0.78 0.75 0.92 0.01 0.33 0.00
6 1 19747.98 20047.94 1.00 1.00 1.00 1.00 1.00
6 2 19380.08 19984.61 0.70 0.89 0.94 0.45 0.55 0.00
6 3 19286.96 20196.06 0.74 0.88 0.89 0.23 0.39 0.00
6 4 19262.78 20476.45 0.78 0.83 0.93 0.11 0.43 1.00
6 5 19282.17 20800.41 0.78 0.81 0.92 0.07 0.43 0.43
6 6 19257.62 21080.44 0.86 0.83 1.00 0.03 0.51 0.00
6 7 19211.10 21338.49 0.84 0.78 1.00 0.01 0.47
6 8 19629.33 22061.29 0.86 0.00 0.95 0.00 0.53
6 9
6 10 19648.29 22675.55 0.87 0.00 1.00 0.00 0.42
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
Compare tidyLPA solutions:
Model Classes AIC BIC
1 1 21342.450 21434.878
1 2 20746.522 20889.786
1 3 20464.673 20658.772
1 4 20341.265 20586.199
1 5 20245.437 20541.207
1 6 20175.056 20521.661
6 1 19889.997 20190.388
6 2 19515.286 20120.691
6 3 19437.181 20347.598
6 4 19389.596 20605.025
6 5 19400.107 20920.550
6 6 19404.469 21229.925
Best model according to AIC is Model 6 with 4 classes.
Best model according to BIC is Model 6 with 2 classes.
An analytic hierarchy process, based on the fit indices AIC, AWE, BIC, CLC, and KIC (Akogul & Erisoglu, 2017), suggests the best solution is Model 6 with 2 classes.
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxubTIgPC0gbHBhX3J3X2RmICU+JVxuICBzY2FsZSgpICU+JVxuICBlc3RpbWF0ZV9wcm9maWxlcygyLFxuICAgIHZhcmlhbmNlcyA9IFwidmFyeWluZ1wiLFxuICAgIGNvdmFyaWFuY2VzID0gXCJ2YXJ5aW5nXCIsXG4gICAgcGFja2FnZSA9IFwiTXBsdXNBdXRvbWF0aW9uXCJcbiAgKVxuYGBgIn0= -->
```r
m2 <- lpa_rw_df %>%
scale() %>%
estimate_profiles(2,
variances = "varying",
covariances = "varying",
package = "MplusAutomation"
)
m2_plot <- lpa_rw_df %>%
scale() %>%
estimate_profiles(2,
variances = "varying",
covariances = "varying"
#package = "MplusAutomation"
) %>%
plot_profiles(add_line = TRUE, rawdata= FALSE, bw = FALSE)
Profile 1: dysfunctional Profile 2: adaptive
get_estimates(m2)
out <- get_data(m2)
lpa_rw_df$lpa_class <- out$Class
table(
lpa_rw_df$lpa_class
)
1 2
332 419
table(
lpa_rw_df$lpa_class, rw_df$job_qualification
)
driver team_leader team_member
1 57 135 140
2 98 164 157
lpa_rw_df$class <- factor(lpa_rw_df$lpa_class)
summary(lpa_rw_df$class)
1 2
332 419
rw_df$class <- lpa_rw_df$class
m1 <- brm(
bf(scs_ts ~ class),
family = student(),
data = rw_df,
init = 0.1,
backend = "cmdstanr",
adapt_delta = 0.9
)
pp_check(m1)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
summary(m1)
Family: student
Links: mu = identity; sigma = identity; nu = identity
Formula: scs_ts ~ class
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 74.19 0.83 72.60 75.82 1.00 3812 3100
class2 13.07 1.10 11.00 15.21 1.00 4001 3053
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 14.30 0.47 13.38 15.20 1.00 3617 2817
nu 28.45 13.18 11.35 62.77 1.00 3428 2922
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
me <- conditional_effects(
m1, "class"
)
plot(me, points = FALSE)
BFt <- BayesFactor::ttestBF(
rw_df$scs_ts[rw_df$class == 1],
rw_df$scs_ts[rw_df$class == 2],
paired = FALSE
)
effectsize(BFt, type = "d")
Cohen's d | 95% CI
--------------------------
-0.87 | [-1.02, -0.72]
- Estimated using pooled SD.
m1 %>%
emmeans( ~ class) %>%
gather_emmeans_draws() %>%
ggplot(aes(x = class, y = .value)) +
geom_eye() +
stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
# facet_grid(~ wool) +
# theme_light()
labs(x = "LPA Class", y = "SCS Score", title = "Rescue Workers") +
papaja::theme_apa() +
annotate("text", x = 1, y = 78.5, label = "Bayesian Cohen's d = 0.89\n 95% CI [0.73, 1.04]")
```r
names(all_items)
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
### Self-judgment
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxubTIgPC0gYnJtKFxuICBiZihzaiB+IGNsYXNzKSxcbiAgZGF0YSA9IHJ3X2RmLFxuICBmYW1pbHkgPSBzdHVkZW50LFxuICBiYWNrZW5kID0gXCJjbWRzdGFuclwiLFxuICBhZGFwdF9kZWx0YSA9IDAuOTlcbilcbmBgYCJ9 -->
```r
m2 <- brm(
bf(sj ~ class),
data = rw_df,
family = student,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m2)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
summary(m2)
Family: student
Links: mu = identity; sigma = identity; nu = identity
Formula: sj ~ class
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 17.12 0.25 16.63 17.62 1.00 3494 2516
class2 -3.67 0.32 -4.29 -3.06 1.00 3404 2787
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 4.26 0.12 4.03 4.51 1.00 3787 2427
nu 46.52 16.90 21.85 86.76 1.00 3495 2719
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
me <- conditional_effects(
m2, "class"
)
plot(me, points = FALSE)
emmeans(m2, specs = pairwise ~ class)
$emmeans
class emmean lower.HPD upper.HPD
1 17.1 16.6 17.6
2 13.4 13.0 13.9
Point estimate displayed: median
HPD interval probability: 0.95
$contrasts
contrast estimate lower.HPD upper.HPD
class1 - class2 3.66 3.07 4.3
Point estimate displayed: median
HPD interval probability: 0.95
BFt <- BayesFactor::ttestBF(
rw_df$sj[rw_df$class == 1],
rw_df$sj[rw_df$class == 2],
paired = FALSE
)
effectsize(BFt)
Cohen's d | 95% CI
------------------------
0.83 | [0.67, 0.98]
- Estimated using pooled SD.
p2 <- m2 %>%
emmeans( ~ class) %>%
gather_emmeans_draws() %>%
ggplot(aes(x = class, y = .value)) +
geom_eye() +
stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
# facet_grid(~ wool) +
# theme_light()
scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
labs(x = "LPA Class", y = "SCS Self-Judgment", title = "A") +
# papaja::theme_apa() +
annotate("text", x = 1, y = 16, label = "Bayesian Cohen's d = 0.82\n 95% CI [0.67, 0.97]")
Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").
p2
m3 <- brm(
bf(is ~ class),
family = skew_normal(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m3)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
summary(m3)
Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: is ~ class
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 12.37 0.25 11.86 12.87 1.00 2367 2295
class2 -3.45 0.33 -4.11 -2.80 1.00 2399 2641
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 4.05 0.12 3.83 4.30 1.00 2133 2366
alpha 2.19 0.66 0.88 3.46 1.00 1966 1683
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
me <- conditional_effects(
m3, "class"
)
plot(me, points = FALSE)
emmeans(m3, specs = pairwise ~ class)
$emmeans
class emmean lower.HPD upper.HPD
1 12.37 11.84 12.8
2 8.93 8.54 9.3
Point estimate displayed: median
HPD interval probability: 0.95
$contrasts
contrast estimate lower.HPD upper.HPD
class1 - class2 3.44 2.81 4.11
Point estimate displayed: median
HPD interval probability: 0.95
BFt <- BayesFactor::ttestBF(
rw_df$is[rw_df$class == 1],
rw_df$is[rw_df$class == 2],
paired = FALSE
)
effectsize(BFt)
Cohen's d | 95% CI
------------------------
0.94 | [0.79, 1.10]
- Estimated using pooled SD.
p3 <- m3 %>%
emmeans( ~ class) %>%
gather_emmeans_draws() %>%
ggplot(aes(x = class, y = .value)) +
geom_eye() +
stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
# facet_grid(~ wool) +
# theme_light()
scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
labs(x = "LPA Class", y = "SCS Isolation", title = "B") +
# papaja::theme_apa() +
annotate("text", x = 1, y = 11, label = "Bayesian Cohen's d = 0.94\n 95% CI [0.79, 1.10]")
Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").
p3
m4 <- brm(
bf(oi ~ class),
family = skew_normal(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m4)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
summary(m4)
Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: oi ~ class
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 11.10 0.25 10.58 11.57 1.00 1434 1663
class2 -2.68 0.35 -3.35 -1.95 1.01 1372 1365
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 3.63 0.12 3.41 3.87 1.00 1445 1791
alpha 3.24 0.91 1.94 5.48 1.01 1314 1122
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
me <- conditional_effects(
m4, "class"
)
plot(me, points = FALSE)
emmeans(m4, specs = pairwise ~ class)
$emmeans
class emmean lower.HPD upper.HPD
1 11.11 10.60 11.58
2 8.42 8.05 8.79
Point estimate displayed: median
HPD interval probability: 0.95
$contrasts
contrast estimate lower.HPD upper.HPD
class1 - class2 2.69 2.01 3.4
Point estimate displayed: median
HPD interval probability: 0.95
BFt <- BayesFactor::ttestBF(
rw_df$oi[rw_df$class == 1],
rw_df$oi[rw_df$class == 2],
paired = FALSE
)
effectsize(BFt)
Cohen's d | 95% CI
------------------------
0.99 | [0.83, 1.15]
- Estimated using pooled SD.
p4 <- m4 %>%
emmeans( ~ class) %>%
gather_emmeans_draws() %>%
ggplot(aes(x = class, y = .value)) +
geom_eye() +
stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
# facet_grid(~ wool) +
# theme_light()
scale_x_discrete(labels=c('High Resilience', 'Low Resilience')) +
labs(x = "LPA Class", y = "SCS Over-Identification", title = "C") +
# papaja::theme_apa() +
annotate("text", x = 1, y = 9.4, label = "Bayesian Cohen's d = 0.97\n 95% CI [0.82, 1.12]")
Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").
p4
m5 <- brm(
bf(sk ~ class),
family = student(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m5)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
summary(m5)
Family: student
Links: mu = identity; sigma = identity; nu = identity
Formula: sk ~ class
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 13.07 0.24 12.59 13.54 1.00 3482 2388
class2 1.18 0.32 0.55 1.80 1.00 3684 2519
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 4.22 0.11 4.00 4.45 1.00 3537 2868
nu 41.60 16.59 17.38 80.94 1.00 3185 2731
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
me <- conditional_effects(
m5, "class"
)
plot(me, points = FALSE)
emmeans(m5, specs = pairwise ~ class)
$emmeans
class emmean lower.HPD upper.HPD
1 13.1 12.6 13.5
2 14.2 13.9 14.7
Point estimate displayed: median
HPD interval probability: 0.95
$contrasts
contrast estimate lower.HPD upper.HPD
class1 - class2 -1.17 -1.79 -0.545
Point estimate displayed: median
HPD interval probability: 0.95
BFt <- BayesFactor::ttestBF(
rw_df$sk[rw_df$class == 1],
rw_df$sk[rw_df$class == 2],
paired = FALSE
)
effectsize(BFt)
Cohen's d | 95% CI
--------------------------
-0.27 | [-0.41, -0.13]
- Estimated using pooled SD.
p5 <- m5 %>%
emmeans( ~ class) %>%
gather_emmeans_draws() %>%
ggplot(aes(x = class, y = .value)) +
geom_eye() +
stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
# facet_grid(~ wool) +
# theme_light()
scale_x_discrete(labels=c('High Resilience', 'Low Resilience')) +
labs(x = "LPA Class", y = "SCS Self-Kindness", title = "D") +
# papaja::theme_apa() +
annotate("text", x = 1, y = 14.5, label = "Bayesian Cohen's d = 0.28\n 95% CI [0.14, 0.43]")
Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").
p5
m6 <- brm(
bf(ch ~ class),
family = student(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m6)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
summary(m6)
Family: student
Links: mu = identity; sigma = identity; nu = identity
Formula: ch ~ class
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 11.61 0.18 11.27 11.96 1.00 3502 2643
class2 -0.04 0.25 -0.53 0.45 1.00 3016 2397
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 3.28 0.09 3.10 3.47 1.00 2680 2104
nu 42.34 16.50 18.20 80.94 1.00 2954 2163
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
me <- conditional_effects(
m6, "class"
)
plot(me, points = FALSE)
BFt <- BayesFactor::ttestBF(
rw_df$ch[rw_df$class == 1],
rw_df$ch[rw_df$class == 2],
paired = FALSE
)
effectsize(BFt)
Cohen's d | 95% CI
-------------------------
9.51e-03 | [-0.13, 0.15]
- Estimated using pooled SD.
p6 <- m6 %>%
emmeans( ~ class) %>%
gather_emmeans_draws() %>%
ggplot(aes(x = class, y = .value)) +
geom_eye() +
stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
# facet_grid(~ wool) +
# theme_light()
scale_x_discrete(labels=c('High Resilience', 'Low Resilience')) +
labs(x = "LPA Class", y = "SCS Common-Humanity", title = "E") +
# papaja::theme_apa() +
annotate("text", x = 1.5, y = 12.3, label = "Bayesian Cohen's d = 0.00\n 95% CI [-0.14, 0.14]")
Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").
p6
m7 <- brm(
bf(mi ~ class),
family = student(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m7)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
summary(m7)
Family: student
Links: mu = identity; sigma = identity; nu = identity
Formula: mi ~ class
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 12.76 0.17 12.43 13.09 1.00 3277 2803
class2 0.96 0.23 0.51 1.42 1.00 3209 2387
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 3.04 0.09 2.87 3.22 1.00 2578 2822
nu 36.60 15.48 15.14 74.03 1.00 2806 3032
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
emmeans(m7, specs = pairwise ~ class)
$emmeans
class emmean lower.HPD upper.HPD
1 12.8 12.4 13.1
2 13.7 13.4 14.0
Point estimate displayed: median
HPD interval probability: 0.95
$contrasts
contrast estimate lower.HPD upper.HPD
class1 - class2 -0.964 -1.43 -0.52
Point estimate displayed: median
HPD interval probability: 0.95
BFt <- BayesFactor::ttestBF(
rw_df$mi[rw_df$class == 1],
rw_df$mi[rw_df$class == 2],
paired = FALSE
)
effectsize(BFt)
Cohen's d | 95% CI
--------------------------
-0.30 | [-0.44, -0.15]
- Estimated using pooled SD.
me <- conditional_effects(
m7, "class"
)
plot(me, points = FALSE)
p7 <- m7 %>%
emmeans( ~ class) %>%
gather_emmeans_draws() %>%
ggplot(aes(x = class, y = .value)) +
geom_eye() +
stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
# facet_grid(~ wool) +
# theme_light()
scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
labs(x = "LPA Class", y = "SCS Mindfulness", title = "F") +
# papaja::theme_apa() +
annotate("text", x = 1, y = 13.8, label = "Bayesian Cohen's d = 0.32\n 95% CI [0.17, 0.46]")
Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").
p7
fig_scs <- (p2 | p3 | p4) /
(p5 | p6 | p7)
out <- fig_scs + plot_annotation(
title = 'SCS Subscales as a Function of LPA Class',
subtitle = 'Rescue Workers group'
# caption = 'Disclaimer: None of these plots are insightful'
)
ggsave("scs_subscales_lpa.pdf", width = 35, height = 20, units = "cm")
m10 <- brm(
bf(ies_ts ~ class),
family = skew_normal(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m10)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
summary(m10)
Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: ies_ts ~ class
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 21.08 0.72 19.77 22.68 1.00 1705 1915
class2 -4.26 1.05 -6.71 -2.57 1.00 1211 1474
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 13.95 0.40 13.19 14.75 1.00 1229 1470
alpha 10.07 2.34 5.90 15.14 1.00 973 1426
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
p10 <- m10 %>%
emmeans( ~ class) %>%
gather_emmeans_draws() %>%
ggplot(aes(x = class, y = .value)) +
geom_eye() +
stat_summary(aes(group = NA), fun.y = mean, geom = "line") +
# facet_grid(~ wool) +
# theme_light()
scale_x_discrete(labels=c('Low Resilience', 'High Resilience')) +
labs(x = "LPA Class", y = "Impact of Event Scale - Revised (IES-R)") +
# papaja::theme_apa() +
annotate("text", x = 1, y = 17, label = "Bayesian Cohen's d = 1.34\n 95% CI [1.18, 1.50]")
Warning: 'geom_eye' is deprecated.
Use 'stat_eye' instead.
See help("Deprecated") and help("tidybayes-deprecated").
p10
BFt <- BayesFactor::ttestBF(
rw_df$ies_ts[rw_df$class == 1],
rw_df$ies_ts[rw_df$class == 2],
paired = FALSE
)
t is large; approximation invoked.
effectsize(BFt)
Cohen's d | 95% CI
------------------------
1.37 | [1.22, 1.53]
- Estimated using pooled SD.
emmeans(m10 , specs = pairwise ~ class)
$emmeans
class emmean lower.HPD upper.HPD
1 21.0 19.6 22.4
2 16.8 15.4 18.2
Point estimate displayed: median
HPD interval probability: 0.95
$contrasts
contrast estimate lower.HPD upper.HPD
class1 - class2 4.13 2.38 6.41
Point estimate displayed: median
HPD interval probability: 0.95
```r
m11 <- brm(
bf(ptgi_total_score | trunc(lb = 0) ~ class),
family = student(),
data = rw_df,
backend = \cmdstanr\
)
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxucHBfY2hlY2sobTExKVxuYGBgXG5gYGAifQ== -->
```r
```r
pp_check(m11)
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuc3VtbWFyeShtMTEpXG5gYGBcbmBgYCJ9 -->
```r
```r
summary(m11)
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
### Job qualification
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxucndfZGYkam9iX3F1YWxpZmljYXRpb24gPC0gaWZlbHNlKFxuICByd19kZiRqb2JfcXVhbGlmaWNhdGlvbiA9PSBcIm5vbl9yZXNjdWVfd29ya2VyXCIsIFwidGVhbV9tZW1iZXJcIiwgXG4gIHJ3X2RmJGpvYl9xdWFsaWZpY2F0aW9uXG4pIFxuYGBgIn0= -->
```r
rw_df$job_qualification <- ifelse(
rw_df$job_qualification == "non_rescue_worker", "team_member",
rw_df$job_qualification
)
m12 <- brm(
bf(ies_ts ~ job_qualification),
family = skew_normal(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m12)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
summary(m12)
Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: ies_ts ~ job_qualification
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept 17.91 0.68 16.56 19.24 1.00
job_qualificationteam_leader 1.20 0.65 -0.02 2.49 1.00
job_qualificationteam_member 1.35 0.66 0.12 2.68 1.00
Bulk_ESS Tail_ESS
Intercept 1914 2193
job_qualificationteam_leader 2133 1952
job_qualificationteam_member 2215 2001
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 14.47 0.39 13.74 15.28 1.00 1942 2043
alpha 15.62 2.23 11.69 20.44 1.00 2406 2098
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
me <- conditional_effects(
m12, "job_qualification"
)
plot(me, points = FALSE)
emmeans(m12 , specs = pairwise ~ job_qualification)
$emmeans
job_qualification emmean lower.HPD upper.HPD
driver 17.9 16.5 19.2
team_leader 19.1 18.0 20.3
team_member 19.2 18.1 20.4
Point estimate displayed: median
HPD interval probability: 0.95
$contrasts
contrast estimate lower.HPD upper.HPD
driver - team_leader -1.182 -2.45 0.0631
driver - team_member -1.345 -2.67 -0.1136
team_leader - team_member -0.149 -1.27 0.8995
Point estimate displayed: median
HPD interval probability: 0.95
m13 <- brm(
bf(scs_ts ~ job_qualification * class),
family = gaussian(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m13)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
me <- conditional_effects(
m13, "job_qualification:class"
)
plot(me, points = FALSE)
summary(m13)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: scs_ts ~ job_qualification * class
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI
Intercept 78.77 1.90 75.10 82.59
job_qualificationteam_leader -3.78 2.30 -8.45 0.61
job_qualificationteam_member -7.34 2.26 -11.84 -3.04
class2 10.15 2.42 5.40 14.80
job_qualificationteam_leader:class2 2.03 2.99 -3.76 7.96
job_qualificationteam_member:class2 4.48 2.94 -1.42 10.19
Rhat Bulk_ESS Tail_ESS
Intercept 1.00 1156 1812
job_qualificationteam_leader 1.00 1185 1852
job_qualificationteam_member 1.00 1235 1860
class2 1.00 1080 1717
job_qualificationteam_leader:class2 1.00 1136 1696
job_qualificationteam_member:class2 1.00 1156 1799
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 14.80 0.39 14.06 15.59 1.00 2650 2594
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
emmeans(m13 , specs = pairwise ~ job_qualification*class)
$emmeans
job_qualification class emmean lower.HPD upper.HPD
driver 1 78.7 75.2 82.6
team_leader 1 75.0 72.4 77.3
team_member 1 71.4 69.0 73.8
driver 2 88.9 86.2 91.9
team_leader 2 87.2 84.9 89.4
team_member 2 86.1 83.8 88.3
Point estimate displayed: median
HPD interval probability: 0.95
$contrasts
contrast estimate lower.HPD upper.HPD
driver class1 - team_leader class1 3.72 -0.306 8.71
driver class1 - team_member class1 7.27 3.137 11.88
driver class1 - driver class2 -10.24 -14.752 -5.37
driver class1 - team_leader class2 -8.43 -12.880 -4.12
driver class1 - team_member class2 -7.34 -11.794 -3.05
team_leader class1 - team_member class1 3.58 0.088 7.06
team_leader class1 - driver class2 -13.93 -18.127 -10.51
team_leader class1 - team_leader class2 -12.17 -15.285 -8.69
team_leader class1 - team_member class2 -11.06 -14.490 -7.82
team_member class1 - driver class2 -17.50 -21.106 -13.43
team_member class1 - team_leader class2 -15.79 -19.155 -12.51
team_member class1 - team_member class2 -14.67 -17.974 -11.21
driver class2 - team_leader class2 1.75 -1.786 5.66
driver class2 - team_member class2 2.85 -1.065 6.51
team_leader class2 - team_member class2 1.12 -1.997 4.29
Point estimate displayed: median
HPD interval probability: 0.95
m14 <- brm(
bf(ptgi_total_score ~ ies_total_score | (1 | comme\)),
family = gaussian(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.99
)
summary(m14)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: ptgi_total_score ~ ies_total_score
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 33.36 1.30 30.82 35.89 1.00 2784 2099
ies_total_score 0.40 0.05 0.29 0.50 1.00 2814 2352
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 23.07 0.61 21.89 24.29 1.00 3030 2390
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
bayes_R2(m14)
Estimate Est.Error Q2.5 Q97.5
R2 0.1089372 0.01913069 0.07305561 0.1471529
conditional_effects(m14, "is:class")
mydf <- data.frame(
scs = scale(rw_df$scs_ts),
class = ifelse(rw_df$class == 1, 0.0, 1.0),
ptgi = scale(rw_df$ptgi_total_score),
psc = scale(rw_df$sk + rw_df$ch + rw_df$mi),
nsc = scale(rw_df$sj + rw_df$oi + rw_df$is),
commettee = rw_df$red_cross_commeetee_location
)
mydf <- mydf[complete.cases(mydf), ]
# Impute NAs with the mode.
# mydf$commettee[is.na(mydf$commettee)] <- "Comitato di Groane"
# summary(factor(mydf$commettee))
f1 <- bf(scs ~ class, family = skew_normal())
f2 <- bf(ptgi ~ scs + class, family = skew_normal())
mod <- brm(
f1 + f2 + set_rescor(FALSE),
data = mydf,
cores = 4,
refresh = 0,
init = 0.1,
backend = "cmdstanr",
adapt_delta = 0.99
)
bayestestR::mediation(mod, mediator = "scs", ci = 0.95, method = "SPI")
# Causal Mediation Analysis for Stan Model
Treatment: class
Mediator : scs
Response : ptgi
Effect | Estimate | 95% SPI
----------------------------------------------------
Direct Effect (ADE) | -0.303 | [-0.461, -0.144]
Indirect Effect (ACME) | 0.120 | [ 0.055, 0.184]
Mediator Effect | 0.153 | [ 0.070, 0.228]
Total Effect | -0.183 | [-0.337, -0.049]
Proportion mediated: -65.54% [-257.06%, 125.97%]
Direct and indirect effects have opposite directions. The proportion mediated is not meaningful.
pp_check(mod, resp = "ptgi")
Using 10 posterior draws for ppc type 'dens_overlay' by default.
pp_check(mod, resp = "scs")
Using 10 posterior draws for ppc type 'dens_overlay' by default.
summary(mod)
Family: MV(skew_normal, skew_normal)
Links: mu = identity; sigma = identity; alpha = identity
mu = identity; sigma = identity; alpha = identity
Formula: scs ~ class
ptgi ~ scs + class
Data: mydf (Number of observations: 746)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
scs_Intercept -0.44 0.05 -0.54 -0.35 1.00 3911 2852
ptgi_Intercept 0.17 0.06 0.05 0.29 1.00 3619 3231
scs_class 0.79 0.07 0.65 0.92 1.00 3666 2804
ptgi_scs 0.15 0.04 0.07 0.23 1.00 3199 2879
ptgi_class -0.30 0.08 -0.47 -0.15 1.00 3176 2768
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_scs 0.91 0.02 0.87 0.96 1.00 4070 2612
sigma_ptgi 0.99 0.03 0.94 1.04 1.00 4207 2703
alpha_scs -0.70 0.59 -1.52 0.59 1.00 2244 3333
alpha_ptgi 0.09 0.53 -0.88 1.10 1.00 3801 3079
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
fit <- mlm(d = mydf,
id = "commettee",
x = "class",
m = "scs",
y = "ptgi",
iter = 2000,
cores = 4)
mlm_path_plot(fit, level = .95, text = T,
xlab = "Resilience\nProfile",
mlab = "Self\nCompassion",
ylab = "PTG", digits = 2)
Treatment: class Mediator : scs Response : ptgi
Direct Effect (ADE) | -0.300 | [-0.464, -0.140] Indirect Effect (ACME) | 0.121 | [ 0.053, 0.187] Mediator Effect | 0.154 | [ 0.069, 0.230] Total Effect | -0.178 | [-0.325, -0.027]
m16 <- brm(
bf(ptgi_total_score ~ class),
family = skew_normal(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m16)
summary(m16)
Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: ptgi_total_score ~ class
Data: rw_df (Number of observations: 751)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 43.19 1.33 40.51 45.70 1.00 3225 2321
class2 -4.36 1.76 -7.72 -0.76 1.00 3158 2300
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 23.85 0.63 22.65 25.19 1.00 3179 2416
alpha 0.39 0.65 -0.73 1.68 1.00 2880 2317
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
bayes_R2(m16)
Estimate Est.Error Q2.5 Q97.5
R2 0.009530306 0.006634892 0.0003035324 0.02528022
m17 <- brm(
bf(ptgi_total_score ~ is_rescue_worker),
family = skew_normal(),
data = all_items,
backend = "cmdstanr",
adapt_delta = 0.99
)
summary(m17)
Family: skew_normal
Links: mu = identity; sigma = identity; alpha = identity
Formula: ptgi_total_score ~ is_rescue_worker
Data: all_items (Number of observations: 1068)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 50.37 1.33 47.76 52.97 1.00 3038 2433
is_rescue_workerSi -9.65 1.60 -12.78 -6.50 1.00 2988 2273
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 23.46 0.52 22.48 24.55 1.00 3236 2596
alpha 0.08 0.47 -0.76 0.98 1.00 3008 2787
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
m19 <- brm(
bf(rate_of_activity_num ~ class),
family = skew_normal(),
data = rw_df,
backend = "cmdstanr",
adapt_delta = 0.99
)
pp_check(m19)
Using 10 posterior draws for ppc type 'dens_overlay' by default.